{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "共享单车数据——线性回归分析\n",
    "美国 Washington,  D.C.的一个共享单车公司提供的自行车数据上进行回归分析，涵盖了16列参数信息。\n",
    "本项目将原始数据集存为csv格式，方便调用pandas做数据分析。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "导入必要的工具包"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np  # 矩阵操作\n",
    "import pandas as pd # SQL数据处理\n",
    "\n",
    "from sklearn.metrics import r2_score  #评价回归预测模型的性能\n",
    "\n",
    "import matplotlib.pyplot as plt   #画图\n",
    "import seaborn as sns\n",
    "\n",
    "# 图形出现在Notebook里而不是新窗口\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "数据探索¶\n",
    "读取数据"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>instant</th>\n",
       "      <th>season</th>\n",
       "      <th>yr</th>\n",
       "      <th>mnth</th>\n",
       "      <th>holiday</th>\n",
       "      <th>weekday</th>\n",
       "      <th>workingday</th>\n",
       "      <th>weathersit</th>\n",
       "      <th>temp</th>\n",
       "      <th>atemp</th>\n",
       "      <th>hum</th>\n",
       "      <th>windspeed</th>\n",
       "      <th>casual</th>\n",
       "      <th>registered</th>\n",
       "      <th>cnt</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>6</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>0.344167</td>\n",
       "      <td>0.363625</td>\n",
       "      <td>0.805833</td>\n",
       "      <td>0.160446</td>\n",
       "      <td>331</td>\n",
       "      <td>654</td>\n",
       "      <td>985</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>0.363478</td>\n",
       "      <td>0.353739</td>\n",
       "      <td>0.696087</td>\n",
       "      <td>0.248539</td>\n",
       "      <td>131</td>\n",
       "      <td>670</td>\n",
       "      <td>801</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0.196364</td>\n",
       "      <td>0.189405</td>\n",
       "      <td>0.437273</td>\n",
       "      <td>0.248309</td>\n",
       "      <td>120</td>\n",
       "      <td>1229</td>\n",
       "      <td>1349</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0.200000</td>\n",
       "      <td>0.212122</td>\n",
       "      <td>0.590435</td>\n",
       "      <td>0.160296</td>\n",
       "      <td>108</td>\n",
       "      <td>1454</td>\n",
       "      <td>1562</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0.226957</td>\n",
       "      <td>0.229270</td>\n",
       "      <td>0.436957</td>\n",
       "      <td>0.186900</td>\n",
       "      <td>82</td>\n",
       "      <td>1518</td>\n",
       "      <td>1600</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   instant  season  yr  mnth  holiday  weekday  workingday  weathersit  \\\n",
       "0        1       1   0     1        0        6           0           2   \n",
       "1        2       1   0     1        0        0           0           2   \n",
       "2        3       1   0     1        0        1           1           1   \n",
       "3        4       1   0     1        0        2           1           1   \n",
       "4        5       1   0     1        0        3           1           1   \n",
       "\n",
       "       temp     atemp       hum  windspeed  casual  registered   cnt  \n",
       "0  0.344167  0.363625  0.805833   0.160446     331         654   985  \n",
       "1  0.363478  0.353739  0.696087   0.248539     131         670   801  \n",
       "2  0.196364  0.189405  0.437273   0.248309     120        1229  1349  \n",
       "3  0.200000  0.212122  0.590435   0.160296     108        1454  1562  \n",
       "4  0.226957  0.229270  0.436957   0.186900      82        1518  1600  "
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# path to where the data lies\n",
    "#dpath = './data/'\n",
    "data = pd.read_csv(\"day.csv\")\n",
    "\n",
    "#通过观察前5行，了解数据每列（特征）的概况\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "数据基本信息\n",
    "样本数目、特征维数 每个特征的类型、空值样本的数目、数据类型"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(731, 15)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "数据探索\n",
    "数据准备"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 从原始数据中分离输入特征x和输出y\n",
    "y = data['cnt'].values\n",
    "X = data.drop('cnt', axis = 1)\n",
    "\n",
    "#用于后续显示权重系数对应的特征\n",
    "columns = X.columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "作业1：训练数据和测试数据分割\n",
    "当数据量比较大时，可用train_test_split从训练集中分出一部分做校验集； 样本数目较少时，建议用交叉验证 在线性回归中，留一交叉验证有简便计算方式，无需显式交叉验证"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(365, 14)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#将数据分割训练数据与测试数据\n",
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "# 随机采样20%的数据构建测试样本，其余作为训练样本\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=50, test_size=0.5)\n",
    "X_train.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "作业2：数据预处理／特征工程\n",
    "特征工程是实际任务中特别重要的环节。\n",
    "scikit learn中提供的数据预处理功能： http://scikit-learn.org/stable/modules/preprocessing.html http://scikit-learn.org/stable/modules/classes.html#module- sklearn.feature_extraction\n",
    "#发现各特征差异较大，需要进行数据标准化预处理\n",
    "#标准化的目的在于避免原始特征值差异过大，导致训练得到的参数权重不归一，无法比较各特征的重要性"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "D:\\Anaconda2\\lib\\site-packages\\sklearn\\utils\\validation.py:475: DataConversionWarning: Data with input dtype int64 was converted to float64 by StandardScaler.\n",
      "  warnings.warn(msg, DataConversionWarning)\n",
      "D:\\Anaconda2\\lib\\site-packages\\sklearn\\utils\\validation.py:475: DataConversionWarning: Data with input dtype int64 was converted to float64 by StandardScaler.\n",
      "  warnings.warn(msg, DataConversionWarning)\n",
      "D:\\Anaconda2\\lib\\site-packages\\sklearn\\utils\\validation.py:475: DataConversionWarning: Data with input dtype int64 was converted to float64 by StandardScaler.\n",
      "  warnings.warn(msg, DataConversionWarning)\n"
     ]
    }
   ],
   "source": [
    "# 数据标准化\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "\n",
    "# 分别初始化对特征和目标值的标准化器\n",
    "ss_X = StandardScaler()\n",
    "ss_y = StandardScaler()\n",
    "\n",
    "# 分别对训练和测试数据的特征以及目标值进行标准化处理\n",
    "X_train = ss_X.fit_transform(X_train)\n",
    "X_test = ss_X.transform(X_test)\n",
    "\n",
    "#对y做标准化不是必须\n",
    "#对y标准化的好处是不同问题的w差异不太大，同时正则参数的范围也有限\n",
    "y_train = ss_y.fit_transform(y_train.reshape(-1, 1))\n",
    "y_test = ss_y.transform(y_test.reshape(-1, 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "尝试缺省参数的线性回归"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>coef</th>\n",
       "      <th>columns</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>[0.7992015607529477]</td>\n",
       "      <td>registered</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>[0.35431900383432696]</td>\n",
       "      <td>casual</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>[2.664146216407351e-15]</td>\n",
       "      <td>instant</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>[4.306471171686414e-16]</td>\n",
       "      <td>temp</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>[9.683329465026209e-17]</td>\n",
       "      <td>hum</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>[6.768049791686516e-18]</td>\n",
       "      <td>weekday</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>[-1.3975909410059654e-17]</td>\n",
       "      <td>weathersit</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>[-8.086278779127036e-17]</td>\n",
       "      <td>windspeed</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>[-1.9648465161928974e-16]</td>\n",
       "      <td>holiday</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>[-2.6278332966459976e-16]</td>\n",
       "      <td>season</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>[-3.1207395989199734e-16]</td>\n",
       "      <td>workingday</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>[-3.7441501753193594e-16]</td>\n",
       "      <td>atemp</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>[-7.378407600702821e-16]</td>\n",
       "      <td>mnth</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>[-2.2732219189246263e-15]</td>\n",
       "      <td>yr</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                         coef     columns\n",
       "13       [0.7992015607529477]  registered\n",
       "12      [0.35431900383432696]      casual\n",
       "0     [2.664146216407351e-15]     instant\n",
       "8     [4.306471171686414e-16]        temp\n",
       "10    [9.683329465026209e-17]         hum\n",
       "5     [6.768049791686516e-18]     weekday\n",
       "7   [-1.3975909410059654e-17]  weathersit\n",
       "11   [-8.086278779127036e-17]   windspeed\n",
       "4   [-1.9648465161928974e-16]     holiday\n",
       "1   [-2.6278332966459976e-16]      season\n",
       "6   [-3.1207395989199734e-16]  workingday\n",
       "9   [-3.7441501753193594e-16]       atemp\n",
       "3    [-7.378407600702821e-16]        mnth\n",
       "2   [-2.2732219189246263e-15]          yr"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 线性回归\n",
    "#class sklearn.linear_model.LinearRegression(fit_intercept=True, normalize=False, copy_X=True, n_jobs=1)\n",
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "# 使用默认配置初始化\n",
    "lr = LinearRegression()\n",
    "\n",
    "# 训练模型参数\n",
    "lr.fit(X_train, y_train)\n",
    "\n",
    "# 预测\n",
    "y_test_pred_lr = lr.predict(X_test)\n",
    "y_train_pred_lr = lr.predict(X_train)\n",
    "\n",
    "\n",
    "# 看看各特征的权重系数，系数的绝对值大小可视为该特征的重要性\n",
    "fs = pd.DataFrame({\"columns\":list(columns), \"coef\":list((lr.coef_.T))})\n",
    "fs.sort_values(by=['coef'],ascending=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "模型评价"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The r2 score of LinearRegression on test is 1.0\n",
      "The r2 score of LinearRegression on train is 1.0\n"
     ]
    }
   ],
   "source": [
    "# 使用r2_score评价模型在测试集和训练集上的性能，并输出评估结果\n",
    "#测试集\n",
    "print 'The r2 score of LinearRegression on test is', r2_score(y_test, y_test_pred_lr)\n",
    "#训练集\n",
    "print 'The r2 score of LinearRegression on train is', r2_score(y_train, y_train_pred_lr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAekAAAF5CAYAAAC7qNLUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAHk9JREFUeJzt3Xu4VXW97/H3V1hCWxRDQVEy0DTvgC0JMo2i1B0m7rMr7XjBR9TS6nHX6aJ1ypW7Tnn01MmjdTa7i5iWqCfTXdudHpRtelIDwysW5GUHEiDlBc0U/J4/5li0gLVYc92Yv8V6v55nPXOOMccc8/sbg8Vn/sYY6zciM5EkSeXZrtEFSJKk9hnSkiQVypCWJKlQhrQkSYUypCVJKpQhLUlSoQxpDXgR8UhETG10HY0UEX8XEb+PiLURMXErfu7aiNi7g9dOj4i7eulznoyId/fGuqStyZDWNq29/5w3/c8/Mw/KzPmdrGdsRGREDO6jUhvtUuBjmTksM3+96YtV21+sQnV5RHw9Igb19EOrz3u8p+uRtlWGtFSAAsL/jcAjnSwzPjOHAe8ATgTO6POqpAHOkNaA17a3HRGTImJBRDwfESsj4uvVYndWj89WvckpEbFdRPzXiHgqIlZFxFURMbzNek+rXlsTEV/Y5HNaIuKGiLg6Ip4HTq8++5cR8WxErIiIyyNi+zbry4g4NyKWRMQLEfGPEbFP9Z7nI+K6tstv0sZ2a42IIRGxFhgEPBARv+tse2XmUuBuYEKb9Q+PiO9WdS+PiC+39rQj4k0R8e8R8VxEPBMRczdp05uq57tExM1VW+4D9mmz3GZHMiJifkScWT3fJyJur7b1MxFxTUTs3MG26GgfS8UxpKWNfRP4ZmbuRC0krqvmH1U97lwdov0lcHr1805gb2AYcDlARBwIfAs4GRgNDAf23OSzZgA3ADsD1wDrgU8AuwJTgGnAuZu851jgLcBk4DPA7Ooz3gAcDHyog3a1W2tm/qXqHUOtp7xP+2//q4jYHzgSWNpm9hxgHfAmYCJwNHBm9do/ArcCrwfGAP+rg1VfAbxMbXudQdd66gF8FdgDOIDa9mjpYNmO9rFUHENaA8FPqt7psxHxLLXw7MirwJsiYtfMXJuZ92xh2ZOBr2fm45m5FrgAOKnq7b0f+JfMvCszXwG+CGw6UP4vM/MnmflaZv45Mxdm5j2ZuS4znwT+idqh5bYuzsznM/MR4GHg1urznwNuoRaQXa21XvdHxIvAYmA+1XaMiN2AvwX+ITNfzMxVwDeAk6r3vUrtcPoemflyZm52MVjV6/574IvVOh6mFvx1ycylmXlb9aVjNfB1Nt92rbqyj6WGMqQ1EJyQmTu3/rB577StWcB+wGMR8auIOG4Ly+4BPNVm+ilgMLBb9drvW1/IzJeANZu8//dtJyJiv4j4aUT8oToE/t+o9arbWtnm+Z/bmR5G+7ZUa70Oq9Z/IvBWYIdq/huBJmBFmy9C/wSMql7/DLWe7n1Ru5K+vR7yyKqettvkqXaWa1dEjIqIa6tD7c8DV7P5tmvVlX0sNZQhLbWRmUsy80PUAuZi4IaI2IHNe8EAT1MLqFZ7UTvkuxJYQe3QLgAR8Tpgl00/bpPpbwOPAftWh2I/Ry3cesOWaq1b1lwH/JLa0QGoBetfgF3bfBnaKTMPqt7zh8w8KzP3AD4MfKv1PHQbq6t63rBJja1erB7/ps283ds8/yq17Xlote1OoYNtt4V9LBXHkJbaiIhTImJkZr4GPFvNXk8tRF6jdj631Y+AT0TEuIgYRq3nOzcz11E71/y+iHhbdTHXl+g8cHcEngfWVud9z+m1hm251u74GnB2ROyemSuonXP+HxGxU3WR2j4R8Q6AiPhARLR+YfkTtTBd33Zlmbke+DHQEhF/U53Tn9nm9dXAcuCUiBhU9cbbnj/fEVhL7cK+PYFPd1T4FvaxVBxDWtrYscAj1RXP3wROqs6jvgR8Bbi7OqQ7Gfge8ANqV34/Qe2ip48DVOeMPw5cS61X/QKwilqPsyOfAv5ztew/A3O3sGxXdVhrd2TmQ8C/89cwPA3YHniUWhDfQO0CMIDDgXurbXozcF5mPtHOaj9G7XD6H4Arge9v8vpZ1eetAQ4C/l+b175E7XD8c8DPqAV+R9rdx1tusdQYkdneUTxJvanqvT5L7VB2ewElSZuxJy31kYh4X3XodgdqI3o9BDzZ2Kok9SeGtNR3ZlC7YOtpYF9qh1U9dCWpbh7uliSpUPakJUkq1FYd1H/XXXfNsWPHbs2PlCSpKAsXLnwmM0fWs+xWDemxY8eyYMGCrfmRkiQVJSLqHk3Pw92SJBXKkJYkqVCGtCRJhdqq56QlSe179dVXWbZsGS+/7Ail24qhQ4cyZswYmpqaur0OQ1qSCrBs2TJ23HFHxo4dS0Rv3fxMjZKZrFmzhmXLljFu3Lhur8fD3ZJUgJdffplddtnFgN5GRAS77LJLj4+MGNKSVAgDetvSG/vTkJYkqVCek5akArW0bP31DRo0iEMOOYR169Yxbtw4fvCDH7Dzzjt3+bPOPPNMPvnJT3LggQduNP/KK69kwYIFXH755V1eJ8CwYcNYu3ZtXctOnTqVSy+9lObm5g3zFixYwFVXXcVll13Wrc9vBHvSkiQAXve617Fo0SIefvhhRowYwRVXXNGt9XznO9/ZLKBL0Nzc3OcBvX79+l5dnyEtSdrMlClTWL58+YbpSy65hMMPP5xDDz2UCy+8EIAXX3yR6dOnM378eA4++GDmzp0L1HqxrUNAf//732e//fbjHe94B3ffffeG9Z1++unccMMNG6aHDRsGwNq1a5k2bRqHHXYYhxxyCDfddNNmta1YsYKjjjqKCRMmcPDBB/OLX/yirjbNnz+f4447DoCWlhbOOOMMpk6dyt57771ReF999dVMmjSJCRMm8OEPf3hD8J5zzjk0Nzdz0EEHbdgGUBvy+qKLLuLtb387119/fV211MvD3ZKkjaxfv5558+Yxa9YsAG699VaWLFnCfffdR2Zy/PHHc+edd7J69Wr22GMPfvaznwHw3HPPbbSeFStWcOGFF7Jw4UKGDx/OO9/5TiZOnLjFzx46dCg33ngjO+20E8888wyTJ0/m+OOP3+girB/+8Iccc8wxfP7zn2f9+vW89NJL3WrnY489xh133MELL7zAm9/8Zs455xyWLl3K3Llzufvuu2lqauLcc8/lmmuu4bTTTuMrX/kKI0aMYP369UybNo0HH3yQQw89dEPdd911V7fq2BJDWpIEwJ///GcmTJjAk08+yVve8hbe8573ALWQvvXWWzcE7Nq1a1myZAlHHnkkn/rUp/jsZz/Lcccdx5FHHrnR+u69916mTp3KyJG1Gz6deOKJ/Pa3v91iDZnJ5z73Oe6880622247li9fzsqVK9l99903LHP44Ydzxhln8Oqrr3LCCScwYcKEbrV3+vTpDBkyhCFDhjBq1ChWrlzJvHnzWLhwIYcffviGbTJq1CgArrvuOmbPns26detYsWIFjz766IaQPvHEE7tVQ2c83C1JAv56Tvqpp57ilVde2XBOOjO54IILWLRoEYsWLWLp0qXMmjWL/fbbj4ULF3LIIYdwwQUXcNFFF222zo7+DGnw4MG89tprG9b/yiuvAHDNNdewevVqFi5cyKJFi9htt902+1vjo446ijvvvJM999yTU089lauuuqpb7R0yZMiG54MGDWLdunVkJjNnztzQ1t/85je0tLTwxBNPcOmllzJv3jwefPBBpk+fvlFdO+ywQ7dq6Iw9aUnFqffK5t6+Alo1w4cP57LLLmPGjBmcc845HHPMMXzhC1/g5JNPZtiwYSxfvpympibWrVvHiBEjOOWUUxg2bBhXXnnlRut561vfynnnnceaNWvYaaeduP766xk/fjxQO4+7cOFCPvjBD3LTTTfx6quvArVD5qNGjaKpqYk77riDp57a/K6OTz31FHvuuSdnnXUWL774Ivfffz+nnXZar7R92rRpzJgxg0984hOMGjWKP/7xj7zwwgs8//zz7LDDDgwfPpyVK1dyyy23MHXq1F75zC0xpCWpQI3+AjJx4kTGjx/Ptddey6mnnsrixYuZMmUKULvI6+qrr2bp0qV8+tOfZrvttqOpqYlvf/vbG61j9OjRtLS0MGXKFEaPHs1hhx224SKss846ixkzZjBp0iSmTZu2oSd68skn8773vY/m5mYmTJjA/vvvv1lt8+fP55JLLqGpqYlhw4Z12JOePn36hnGzp0yZwkc/+tFO233ggQfy5S9/maOPPprXXnuNpqYmrrjiCiZPnszEiRM56KCD2HvvvTniiCPq35g9EJm5VT4IoLm5OVuv+JOkjgzEnvTixYs54IADGl2Gell7+zUiFmZmcwdv2YjnpCVJKpQhLUlSoQxpSSrE1jz9qL7XG/vTkJakAgwdOpQ1a9YY1NuI1vtJDx06tEfr8epuSSrAmDFjWLZsGatXr250KeolQ4cOZcyYMT1ahyEtSQVoampi3LhxjS5DhfFwtyRJhTKkJUkqlCEtSVKhDGlJkgplSEuSVChDWpKkQhnSkiQVypCWJKlQhrQkSYUypCVJKlRdw4JGxJPAC8B6YF1mNkfECGAuMBZ4EvhgZv6pb8qUJGng6UpP+p2ZOSEzm6vp84F5mbkvMK+aliRJvaQnh7tnAHOq53OAE3pejiRJalVvSCdwa0QsjIizq3m7ZeYKgOpxVHtvjIizI2JBRCzwFmySJNWv3ltVHpGZT0fEKOC2iHis3g/IzNnAbIDm5mbvZi5JUp3q6kln5tPV4yrgRmASsDIiRgNUj6v6qkhJkgaiTkM6InaIiB1bnwNHAw8DNwMzq8VmAjf1VZGSJA1E9Rzu3g24MSJal/9hZv5bRPwKuC4iZgH/AXyg78qUJGng6TSkM/NxYHw789cA0/qiKEmS5IhjkiQVy5CWJKlQhrQkSYUypCVJKpQhLUlSoQxpSZIKZUhLklQoQ1qSpEIZ0pIkFcqQliSpUIa0JEmFMqQlSSqUIS1JUqEMaUmSCmVIS5JUKENakqRCGdKSJBXKkJYkqVCGtCRJhTKkJUkqlCEtSVKhDGlJkgplSEuSVChDWpKkQhnSkiQVypCWJKlQhrQkSYUypCVJKpQhLUlSoQxpSZIKZUhLklQoQ1qSpEIZ0pIkFcqQliSpUIa0JEmFMqQlSSqUIS1JUqEMaUmSCmVIS5JUqMGNLkCS+qOWlt5dTmqPPWlJkgplSEuSVChDWpKkQhnSkiQVypCWJKlQhrQkSYUypCVJKlTdIR0RgyLi1xHx02p6XETcGxFLImJuRGzfd2VKkjTwdKUnfR6wuM30xcA3MnNf4E/ArN4sTJKkga6ukI6IMcB04DvVdADvAm6oFpkDnNAXBUqSNFDVOyzo/wQ+A+xYTe8CPJuZ66rpZcCe7b0xIs4GzgbYa6+9ul+ppCI5PKbUdzrtSUfEccCqzFzYdnY7i2Z778/M2ZnZnJnNI0eO7GaZkiQNPPX0pI8Ajo+I9wJDgZ2o9ax3jojBVW96DPB035UpSdLA02lPOjMvyMwxmTkWOAm4PTNPBu4A3l8tNhO4qc+qlCRpAOrJ30l/FvhkRCyldo76u71TkiRJgi7eTzoz5wPzq+ePA5N6vyRJkgSOOCZJUrEMaUmSCmVIS5JUKENakqRCGdKSJBWqS1d3S1J/5NCl6q/sSUuSVChDWpKkQhnSkiQVypCWJKlQhrQkSYUypCVJKpQhLUlSoQxpSZIKZUhLklQoQ1qSpEIZ0pIkFcqQliSpUIa0JEmFMqQlSSqUIS1JUqEMaUmSCmVIS5JUKENakqRCGdKSJBXKkJYkqVCGtCRJhTKkJUkqlCEtSVKhDGlJkgplSEuSVChDWpKkQhnSkiQVanCjC5C0dbW09O5ykvqOPWlJkgplSEuSVChDWpKkQhnSkiQVypCWJKlQhrQkSYUypCVJKpQhLUlSoQxpSZIKZUhLklQoQ1qSpEIZ0pIkFcqQliSpUJ2GdEQMjYj7IuKBiHgkIr5UzR8XEfdGxJKImBsR2/d9uZIkDRz19KT/ArwrM8cDE4BjI2IycDHwjczcF/gTMKvvypQkaeDpNKSzZm012VT9JPAu4IZq/hzghD6pUJKkAaquc9IRMSgiFgGrgNuA3wHPZua6apFlwJ4dvPfsiFgQEQtWr17dGzVLkjQg1BXSmbk+MycAY4BJwAHtLdbBe2dnZnNmNo8cObL7lUqSNMB06eruzHwWmA9MBnaOiMHVS2OAp3u3NEmSBrZ6ru4eGRE7V89fB7wbWAzcAby/WmwmcFNfFSlJ0kA0uPNFGA3MiYhB1EL9usz8aUQ8ClwbEV8Gfg18tw/rlCRpwOk0pDPzQWBiO/Mfp3Z+WpIk9QFHHJMkqVCGtCRJhTKkJUkqlCEtSVKhDGlJkgplSEuSVChDWpKkQhnSkiQVypCWJKlQhrQkSYUypCVJKpQhLUlSoQxpSZIKZUhLklQoQ1qSpEIZ0pIkFcqQliSpUIa0JEmFMqQlSSqUIS1JUqEMaUmSCmVIS5JUKENakqRCGdKSJBXKkJYkqVCGtCRJhTKkJUkqlCEtSVKhDGlJkgplSEuSVChDWpKkQhnSkiQVypCWJKlQhrQkSYUypCVJKpQhLUlSoQxpSZIKZUhLklQoQ1qSpEINbnQBktRdLS2NrkDqW/akJUkqlCEtSVKhDGlJkgplSEuSVChDWpKkQhnSkiQVypCWJKlQnYZ0RLwhIu6IiMUR8UhEnFfNHxERt0XEkurx9X1friRJA0c9Pel1wH/JzAOAycBHI+JA4HxgXmbuC8yrpiVJUi/pNKQzc0Vm3l89fwFYDOwJzADmVIvNAU7oqyIlSRqIunROOiLGAhOBe4HdMnMF1IIcGNXBe86OiAURsWD16tU9q1aSpAGk7pCOiGHA/wH+ITOfr/d9mTk7M5szs3nkyJHdqVGSpAGprpCOiCZqAX1NZv64mr0yIkZXr48GVvVNiZIkDUz1XN0dwHeBxZn59TYv3QzMrJ7PBG7q/fIkSRq46rlV5RHAqcBDEbGomvc54GvAdRExC/gP4AN9U6IkSQNTpyGdmXcB0cHL03q3HEmS1MoRxyRJKpQhLUlSoQxpSZIKZUhLklQoQ1qSpELV8ydYkgaglpZGVyDJnrQkSYUypCVJKpQhLUlSoQxpSZIKZUhLklQor+6WtFV4tbjUdfakJUkqlCEtSVKhDGlJkgplSEuSVChDWpKkQhnSkiQVypCWJKlQhrQkSYUypCVJKpQhLUlSoRwWVCpYvUNpOuTmwOK/i4HDnrQkSYUypCVJKpQhLUlSoQxpSZIKZUhLklQoQ1qSpEIZ0pIkFcqQliSpUIa0JEmFMqQlSSqUIS1JUqEMaUmSCmVIS5JUKENakqRCGdKSJBXKkJYkqVCGtCRJhTKkJUkqlCEtSVKhDGlJkgplSEuSVKjBjS5AkkrR0tLoCqSN2ZOWJKlQhrQkSYXqNKQj4nsRsSoiHm4zb0RE3BYRS6rH1/dtmZIkDTz19KSvBI7dZN75wLzM3BeYV01LkqRe1GlIZ+adwB83mT0DmFM9nwOc0Mt1SZI04HX36u7dMnMFQGauiIhRHS0YEWcDZwPstdde3fw4Sdq2eWW52tPnF45l5uzMbM7M5pEjR/b1x0mStM3obkivjIjRANXjqt4rSZIkQfdD+mZgZvV8JnBT75QjSZJa1fMnWD8Cfgm8OSKWRcQs4GvAeyJiCfCealqSJPWiTi8cy8wPdfDStF6uRZIkteGIY5IkFcqQliSpUIa0JEmFMqQlSSqUIS1JUqEMaUmSCmVIS5JUKENakqRCGdKSJBXKkJYkqVCGtCRJhep07G5JUve1tDS6AvVn9qQlSSqUIS1JUqEMaUmSCmVIS5JUKC8ckxrAi4kk1cOetCRJhTKkJUkqlCEtSVKhDGlJkgplSEuSVCiv7pa2AV4trp6o99+P/862PnvSkiQVypCWJKlQhrQkSYUypCVJKpQXjqmhunIhihetSBpo7ElLklQoQ1qSpEIZ0pIkFcqQliSpUIa0JEmF8upubXNXWPf2EIfb2vbRwNHIf48ONdo77ElLklQoQ1qSpEIZ0pIkFcqQliSpUIa0JEmFMqQlSSqUIS1JUqEMaUmSCmVIS5JUKENakqRC9ethQfvDcI39ocb+ore3T19sb/eh1DWNHMa3Xo38vbYnLUlSoQxpSZIK1aOQjohjI+I3EbE0Is7vraIkSVIPQjoiBgFXAH8LHAh8KCIO7K3CJEka6HrSk54ELM3MxzPzFeBaYEbvlCVJkiIzu/fGiPcDx2bmmdX0qcBbM/Njmyx3NnB2Nflm4Dfd+LhdgWe6VWiZtqX2bEttgW2rPbalXNtSe2xL170xM0fWs2BP/gQr2pm3WeJn5mxgdg8+h4hYkJnNPVlHSbal9mxLbYFtqz22pVzbUntsS9/qyeHuZcAb2kyPAZ7uWTmSJKlVT0L6V8C+ETEuIrYHTgJu7p2yJElStw93Z+a6iPgY8HNgEPC9zHyk1yrbWI8OlxdoW2rPttQW2LbaY1vKtS21x7b0oW5fOCZJkvqWI45JklQoQ1qSpEIVGdIRcUlEPBYRD0bEjRGxcwfLPRkRD0XEoohYsLXrrFcX2lP8MKsR8YGIeCQiXouIDv9UoR/tm3rb0x/2zYiIuC0illSPr+9gufXVflkUEUVd7NnZdo6IIRExt3r93ogYu/WrrE8dbTk9Ila32RdnNqLOekTE9yJiVUQ83MHrERGXVW19MCIO29o11quOtkyNiOfa7Jcvbu0aN5KZxf0ARwODq+cXAxd3sNyTwK6Nrrc32kPt4rvfAXsD2wMPAAc2uvZ26jyA2qA084HmLSzXX/ZNp+3pR/vmvwPnV8/P38LvzdpG19rd7QycC/zv6vlJwNxG192DtpwOXN7oWutsz1HAYcDDHbz+XuAWauNnTAbubXTNPWjLVOCnja6z9afInnRm3pqZ66rJe6j9DXa/VWd7+sUwq5m5ODO7M2pckepsT7/YN9RqmlM9nwOc0MBauqOe7dy2jTcA0yKivYGVGq2//JupS2beCfxxC4vMAK7KmnuAnSNi9NaprmvqaEtRigzpTZxB7RtaexK4NSIWVsOP9gcdtWdP4PdtppdV8/qr/rhvOtJf9s1umbkCoHoc1cFyQyNiQUTcExElBXk923nDMtUX3+eAXbZKdV1T77+Zv68OD98QEW9o5/X+or/8jtRrSkQ8EBG3RMRBjSykJ8OC9khE/F9g93Ze+nxm3lQt83lgHXBNB6s5IjOfjohRwG0R8Vj1LWmr64X21DXM6tZQT1vq0K/2TWeraGdecfumC6vZq9o3ewO3R8RDmfm73qmwR+rZzsXsi07UU+e/AD/KzL9ExEeoHSF4V59X1jf6y36px/3UxtZeGxHvBX4C7NuoYhoW0pn57i29HhEzgeOAaVmdKGhnHU9Xj6si4kZqh5gaEgS90J5ihlntrC11rqPf7Js69It9ExErI2J0Zq6oDjWu6mAdrfvm8YiYD0ykdv600erZzq3LLIuIwcBwyjx02WlbMnNNm8l/pna9Sn9VzO9IT2Xm822e/2tEfCsids3MhtxEpMjD3RFxLPBZ4PjMfKmDZXaIiB1bn1O7OKvdq/UarZ72sA0Ns9qf9k2d+su+uRmYWT2fCWx2lCAiXh8RQ6rnuwJHAI9utQq3rJ7t3LaN7wdu7+hLfIN12pZNztkeDyzeivX1tpuB06qrvCcDz7WeeulvImL31uscImIStZxcs+V39aFGX7nW3g+wlNr5jUXVT+vVnHsA/1o935vaFZMPAI9QO3TZ8Nq7255q+r3Ab6n1aopsD/B31L41/wVYCfy8n++bTtvTj/bNLsA8YEn1OKKa3wx8p3r+NuChat88BMxqdN2btGGz7QxcRO0LLsBQ4Prqd+o+YO9G19yDtny1+v14ALgD2L/RNW+hLT8CVgCvVr8vs4CPAB+pXg/giqqtD7GFv/xo9E8dbflYm/1yD/C2RtbrsKCSJBWqyMPdkiTJkJYkqViGtCRJhTKkJUkqlCEtSRpQOrvJRjfW928R8WxE/HST+VdGxBNtbtYxoavrNqQlSQPNlcCxvbi+S4BTO3jt05k5ofpZ1NUVG9KSpAEl27nJRkTsU/WIF0bELyJi/y6sbx7wQm/XCYa0JEkAs4GPZ+ZbgE8B3+ql9X6luonKN1pH++uKho3dLUlSCSJiGLXR+K5vc+fT1uFz/xO1keI2tTwzj+lk1RcAf6B2T/HZ1IaHbm9dHTKkJUkD3XbAs5m52YVdmflj4MfdWWn+dfzyv0TE96n10LtcmCRJA1bW7nz1RER8AKC6Ucj4nq639SYq1Q07TqAbNxpy7G5J0oASET8CpgK7UruxzoXA7cC3gdFAE3BtZtZ1aDoifgHsDwyjdsesWZn584i4HRhJ7QYki6jdxGNtl2o1pCVJKpOHuyVJKpQhLUlSoQxpSZIKZUhLklQoQ1qSpEIZ0pIkFcqQliSpUP8f7ZoNFETQTUEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 504x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#在训练集上观察预测残差的分布，看是否符合模型假设：噪声为0均值的高斯噪声\n",
    "f, ax = plt.subplots(figsize=(7, 5)) \n",
    "f.tight_layout() \n",
    "ax.hist(y_train - y_train_pred_lr,bins=40, label='Residuals Linear', color='b', alpha=.5); \n",
    "ax.set_title(\"Histogram of Residuals\") \n",
    "ax.legend(loc='best');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARgAAADQCAYAAADcQn7hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xt4VNXV+PHvyuQekhBIALkICpSKXOSmeMcLFS+oUBF4qY+0SuWtVi1qvVdtXx+rKFprXyvipVq1imh+AgoCVusN5B5ERAEViQiEACGBJJNk/f6YCe9kZkgGmDNnMlmf58nDnOScs9cTYbn3PvusLaqKMcY4IcntAIwxicsSjDHGMZZgjDGOsQRjjHGMJRhjjGMswRhjHGMJxhjjGEswxhjHWIIxxjgm2e0ADkV+fr5269bN7TCMafGWL19eoqoFTZ3XrBJMt27dWLZsmdthGNPiich3kZxnQyRjjGMswRhjHONaghGRdBH5TERWi8haEbnPrViMMc5wcw6mCjhbVctFJAX4SETeUdXFLsZkjIki13ow6lPuP0zxf1lxGmNcUF1d7ch9XZ2DERGPiKwCtgMLVHVJmHN+LSLLRGTZjh07Yh+kMQmuqKiIQYMG8fLLL0f93q4mGFWtVdUTgM7AiSLSJ8w501V1sKoOLiho8rG7MSZCdXV1TJs2jSFDhlBSUkLbtm2j3kZcPEVS1d3A+8AIl0MxpkXYsmULw4cP56abbuL888+nqKiI8847L+rtuPkUqUBEWvs/ZwDnAl+6FY8xLcnSpUtZsmQJM2bM4M0338Sp0YGbT5GOAv4hIh58ie41VZ3jYjzGJLQ9e/bwySefcP755zNq1Cg2bdpEu3btHG3TtQSjqkXAALfaN6Yl+fDDD7niiivYsWMH3333Hfn5+Y4nF4iTORhjjDOqq6u58847GTZsGB6Ph4ULF5Kfnx+z9pvVy47GmMh5vV5OP/10PvvsM6666ioeffRRsrOzYxqD9WCMSVApKSmMGjWKWbNmMWPGjJgnF7AEY0xC2bZtG5dccgnvvfceALfddhujR492LR5LMMYkiNmzZ9O3b1/mz5/P5s2b3Q4HsARjTLNXUVHB5MmTufjii+nYsSPLly9n4sSJbocFWIIxptl79dVXmT59Or///e9ZsmQJxx9/vNshHWBPkYxphmpra1m3bh19+vRh4sSJ9O/fn0GDBrkdVgjrwRjTzHzzzTeceeaZnH766ZSUlJCUlBSXyQUswRjTbKgqL7zwAv3792fNmjU88cQTjrwBHU02RDKmGfB6vUyYMIGZM2dyxhln8MILL9C1a1e3w2qSJRhj4kzhymKmzl/PD7v307F1Brec14tLB3QiOzubBx54gFtuuQWPx+N2mBGxBGNMHClcWcztb6xhv7cWranm8zefZco3w2HySGbMmIGIuB3iIbEEY0wcmTp/Pfu9tVTv+JaS2Q/j3fEtSZm5TJ3/Ey4d0Mnt8A6ZJRhj4kjxrgrKlr3Frg+eJymtFe0uu4eM7kP4Yfd+t0M7LK4lGBHpArwAdADqgOmq+he34jEm1sLNtXi+/je73ptBRo+TaDvit3iyWgPQsXWGy9EeHjd7MDXATaq6QkSygeUiskBVv3AxJmNionBlMbfMXI23zrdTz+at27hlZiVjxk/gRdJI7nnqgfmWjBQPt5zXy81wD5ub+yJtVdUV/s97gXVA8xtkGnMYbn+jCG+dUle1j5K5j7L1+Ruo2l/OvC928re7r6NzXiYCdGqdwQOj+zbL+ReIkzkYEemGr3xm2H2RgF8DHH300TGNy5houqtwDa8s+Z5a9fVaKrespWTONGrLdpB78lgkOY3d+71cOqBTs00owVxPMCLSCpgF3KiqZcE/V9XpwHSAwYMH286Pplm6q3AN/1zsK6GgdbXs/uhlyhbPJDm3HR0mPEhap+NcjtAZriYY/57Us4CXVPUNN2MxximFK4sPJBcARKj+YR1Zfc6mzTm/Jikt88CP8jJTXIjQOW4+RRLgGWCdqk5zKw5jnFS/cE5VKS96l4xjB5GcnU+7y+5FklNDzr9nZPyUWogGN3swpwJXAGv8+1MD3KGqb7sYkzFRETjfUluxm53zHmf/hs/IOXkseWdcEZJcBJgw9OiEmXup5+a+SB/h+70ak1AC51v2bVzKzrf/Ql1VBXnnTCJ70MiQ8zsFvG+UaFyf5DUmUUx4+lM+3lh64Lh8zSJ2vv0oKQXdaD/uf0gt6NbgfE+S8MiY/gmZWOpZgjHmCAX2WMD3lEiSPGT2PIma0yaQe9JlSHLDydu8zBTuGXl8QicXsARjzBHpd888yqpqAV9iKVsyi31fL6bDhAdJSm9F61PHNzg/kYdD4ViCMeYwBA+HavZso2TOI1Rt+YLMn56O1ngRT8Neyy+GHs3/XNo31qG6yhKMMYeox+1zqfEv+VRVKtb+m9IFTwJC24tuIqv3sAZ1WzwijD+pS4tLLmAJxphDEphcAKiroWzJLFLbHUP+RTeRnNu+wfmndm/DS5NOjm2QccQSjDERKFxZzK2zig4kl8rNa0ht352ktEzajf0TnsxcJKlhGcv22aktOrlABG9Ti0h7EXlGRN7xH/cWkaucD82Y+DDh6U+58dVVVNXUoTVedr33DNteuZ09n74GQHKrNiHJpWe7LJbcOdyNcONKJD2Y54HngDv9x18Br+Jb5m9Mwgp+/BxYxrLVgAvIPWVc2Ota4mTuwUSSYPJV9TURuR1AVWtEpNbhuIxx1fBp7/P19ooDx/vWf8KO2VNJSsui4LJ7yOw+JOSanu2yWDBlWAyjjH+RJJgKEWkLKICIDAX2OBqVMS4J7rXUSz2qB5m9TqHN2ZMOlLEMZMklvEgSzBTgLaC7iHwMFACXORqVMS446f4FbNtbfeC4Yv3H7P/qU9pedBPJOe0oGHlLyDVpyUk8+PN+LWbh3KFqMsH4a+aeCfTC93LielX1Oh6ZMTHiq4+7Cm+d77iuah+lC6dT8flCUjv0pK6yHE9Gdsh1j409wRJLE5pMMCJyLb6CUGv9x3kiMl5V/9fx6Ixx2F2Fa3hp8Wbql7ZUblnHzrmPULNnO7knjyX31PGIp+E/k5QkmDrGkkskIhkiTVLVv9UfqOouEZkEWIIxzVrwcn+t9VIyeyqI0P6//kx6594h17T0hXOHKpIEkyQioqr1k7weILQU12EQkWeBi4DtqtonGvc0JhKBycW7+0eSs/MRTwrtfn43ybntG5SxrGfJ5dBFsm3JfOA1ETlHRM4GXgHmRan954ERUbqXMU0qXFnMgD++y8cbS1FV9q6ax9Znr2XP4pkApLY7JmxyeWzsCZZcDkMkPZhbgWuA/8Y3yfsuMCMajavqf/xblhjjqMKVxdz71lp27/c9n6jdt4ed7zzO/g1LSO96Aq36hV91a4vmjkwkT5HqgCf9X8Y0O8GL5iq/K2LH7Ieoq6wg7+xJZA8eiUjDzvy3f74w1mEmpIMmGBF5TVUvF5E1QMh+RKraz9HI/i8O23jNHLbgtS0ASemtSM4uoO3Y0DKW4JtrMdEh/rnb0B+IHKWqW0Wka7ifq+p3UQnAN0SaE8kk7+DBg3XZsmXRaNYkuODEUvXjBvZvXHqgwpyqNqjZUs8mciMjIstVdXBT5x20B+NPLh7gGVU9N6rRGeOgkDKWn73B7g//iSezNdkDL8STkdMgubTOSOHeixO/Pq4bGp2DUdVaEdknIrmqGvX3j0TkFWAYkC8iW4B7VNXe0jaHpXBlMb9/fTXVtb5eec2e7ZTMnUbV95+T2etU2px3XciKXFuN66xIniJV4tscbQFwYKZMVa8/0sZVdXzTZxnTtHCL5n586Vbqqsppe+HvyDr+7JAh0S8ScKOzeBNJgpnr/zImLgUOieqq9iGpGYgnhbYjriO5TSdSWncIucYeP8dGJI+p/yEiqcBP8T1NWq+q1U1cZkxM/PTOt6n0D4kqNxdRMudRck8eQ/aAC8g4dlDI+TaJG1uRvOx4AfAUsBHfQrtjROQaVX3H6eCMOZjAN6C1xsvuD1+k7LM3SW7TkdQOPcJeY/MtsRfJEGkacJaqbgAQke74hkyWYEzMFa4s5sZXVx04ri7Z7CtjuX0TrU4YQd5ZV5OUmt7gGuu1uCeSBLO9Prn4bQK2OxSPMQcVrtpcbdkOastLKfj53WT2OCnkGksu7ookwawVkbeB1/DNwYwBlorIaABVfcPB+IwJ3UWxvJSqzUVk9R5GxrGD6HTNjJBeC9iQKB5EkmDSgW3Amf7jHUAbYCS+hGMJxjgmOLns++oTds57Aq31kn7MQDwZOSHJJd0jfHn/BbEO1YQRyVOkX8YiEGOCBb6kWFe1j9JFT1OxZgGpHXqQf9HNeDJyQq6xIVF8sZ0dTdwJnsjVWi9bX5hCTWkxOSdfTutTx4dsLJ+T5qHoPistFG8swZi4EviSomodIkmIJ4WcIZeS0rYz6V1C34ltn51quyjGKUswJi4EPyHy7vqBkjmPkHvyWDJ7nEj2CeF7JzaRG98aqwczpbELVXVa9MMxLVGDt59VKS9awK5F0337PdcdfBNRSy7xr7EeTP1rp72AIfg2XwPf06P/OBmUaTm63fZ/r7nV7tvDznl/Zf/Xi0nv2o+2F0whOSc/5BrbRbH5aKwezH0AIvIuMFBV9/qP7wVmxiQ6k7CCy1gC7N+0nP2blpF31q/IHnJpSBlLsF5LcxPJHMzRQODLjdVAN0eiMS1C4AuKdd4qqrdtJL1zb7KOP4u0zr3Dvv1sT4map0gSzIvAZyLyJr6FdaOAFxyNyiSk4Inc6m0bKZn9MDV7S+g0+Rk8GTlhk4v1WpqvSBba3S8i7wCn+7/1S1VdGY3GRWQE8BfAA8xQ1T9H474m/gSuyPWVsXzTX8Yyh4JL7wi7aM5W5DZ/kT6mzgTKVPU5ESkQkWNU9Zsjadhf7/dvwHBgC773m95S1S+O5L4mvgTPtWitl22v/YGqzWvI/MkptBlxna3ITWCR1IO5BxiM72nSc0AK8E/g1CNs+0Rgg6pu8rfzL+ASwBJMggi3ZYh4Ukg76ie06nMOWX3OCVvZ34ZEiSOSHswoYACwAkBVfxCR7MYviUgn4PuA4y1AyPv2ti9S89RgIreynNJF08kecCFpHXuRNyz86202JEo8kSSYalVVEVEAEcmKUtuh/+sKv8HbdGA6+PZFilLbxiHBbz9Xbl5Dydxp1O7dSVqn40jr2CvsdTYkSkyRJJjXROQpoLWITAJ+RXT2pt4CdAk47gz8EIX7Gpc0WJFb62X3hy9RtmQWyXkd6PCLqWGTiwDf2DatCSuSp0gPi8hwoAzfPMwfVHVBFNpeCvQUkWOAYmAc8F9RuK+JseC3nwHKixZQtuR1WvU/j7yzryYpNSPkOuu1JL5IJnkfVNVbgQVhvnfYVLVGRK4D5uN7TP2sqq49knua2At8SqSq1JZtJzm3Pa36n0dym05kdO0f9jqbyG0ZIhkiDQeCk8n5Yb53yFT1beDtI72Pib3gXktNeSk73/4L1ds30vGqJ/FkZIdNLrYfUcvS2NvU/w38BuguIkUBP8oGPnE6MBO/gpPLvq8+Zee8v6LeSvLOuoqk9FZhr/vW5lpanMZ6MC/j25rkAeC2gO/vVdXS8JeYRNegIFStl9J3n6S86F1S23cn/6KbScnvEnKNTeS2XI29Tb0H2CMifwFKA96mzhaRk1R1SayCNO4Lt2UIScnUVu4lZ+hltD5tQkgZS7DSCi1dJHMwTwIDA44rwnzPJLDA5KJ1tZQtfp3M3meS0roDBZfeHrasgiUWA5ElGFHVAwvcVLVORKzUZgsR+JTIu2srJXMepvqH9SBC7smXW80W06hIEsUmEbkeX68FfBO/m5wLycSDwIlcVaVizUJKF00HSSJ/5C1k9T4z7HWWXEygSBLMZOBx4C58S/kX4X83yCSewpXF/O7VVQ3e2ShfPY/S+X8j7ei+5F/4O5Jz2oVcZ0MiE04kK3m341tlaxJc8HtEdd5KklLSyep9FqjS6oQRIUMie0JkGtPYOpjfq+pDIvJXwr+EeL2jkZmYarCLoreK3R88T+W3q+hw5aMkpaaTPSD0LWfbj8g0pbEezDr/n8tiEYhxR+HKYm6ZuQpvne+4etsmSmY/jHfnZrIHXRx2EhdsRa6JTGPrYGb7//xH7MIxsdRwIrfOV8byPy/iycyh3eV/JOOY8CsRbCLXRKqxIdJswgyN6qnqxY5EZGIieL6Fujr2rf+YzB4nHrSMpfVazKFqbIj0sP/P0UAHfGUyAcYD3zoYk3FQ8JCo4suPSO/aD09GDu3H/glJzQwpY2mV5szhamyI9AGAiPxJVc8I+NFsEbGdHZuhwF5LXVUFpe8+ScUX75Nz8ljyzriCpLTQYoVWs8UciUjWwRSIyLEBxbmPAQqcDctE212Faw4kl8rvP6dkziPU7t1J7mkTyD358pDzU5Jg6hibazFHJpIE8zvgfRGpX73bDbjmSBoVkTHAvcBxwImqak+qHBK8bUj52n+zc840klt3oMOEh0jr9NOQa2zRnImWSBbazRORnkD938QvVbXqCNv9HN/czlNHeB9zEME1W1QVESGj2wlkD76Y1qf/wspYGsdFUjIzE5gCdFXVSSLSU0R6qeqcw21UVdf57324tzCNCC5jWb5yLvu+XkK7MffiycqjzTmTQq6xxGKcEMkQ6TlgOVD/t28LMBM47ARzKGxfpEMTWBCqtnwXJe88RuWm5aQfMwj1ViJBE7lJAtMut7kW44xIEkx3VR0rIuMBVHW/RND1EJGF+B5vB7tTVf9fpAHavkiRCV7Xsu/rxex853HUW0mb4ZNpNeDCkB5jVqqH+0f1teRiHBPRxmsikoF/0Z2IdAeanINR1XOPMDYTocD9iAC0xsuu957Bk51P/sibSc1v2POzSVwTK5EkmHuAeUAXEXkJ357UE50MykQmuIxl1Y8bSGnbhaSUNNqN/RPJ2W1DylgKWHIxMdNogvEPhb7E98RnKL6/nzeoasmRNCoio4C/4ltPM1dEVqnqeUdyz5amQfHtulr2LJ7Jno9eJnfoGFqfcQUprUNHp/b2s4m1RhOMf0/qQlUdBMyNVqOq+ibwZrTu15IE91q8u39k55xHqCpeR+ZxZ5Jz4qiQayyxGLdEMkRaLCJDVHWp49GYRoVM5G74jJLZU/1lLG8mq/ewkGssuRg3RZJgzgImi8i3+HYUEHydm35OBmYaCl6RC5CSdxRpnXvT9mfXkpxrZSxN/IkkwZzveBSmUT1un0uN/wH9/m9WsH/jUvLO+TUpbbvQfsx9Ied7BB6xtS0mDjRWDyYdX8HvHsAa4BlVrYlVYKbhkEhrqtn1/vPsXf4WKW27UFdVgSfMFq1Ws8XEk8Z6MP8AvMCH+HoxvYEbYhFUSxf8HlH19m8omT0Vb8lmsgeNpPWZE0lKSWtwjQCPWqU5E2caSzC9VbUvgIg8A3wWm5BatuCJXK3xsn3mvaB1tBtzHxnHDgq5xnotJl41lmC89R9UtcZeTHRe4ERuTXkpnqzWSHIK+ZfcRkqbjngycxucb5O4Jt41lmD6i0iZ/7MAGf7j+qdIoUVbzWEJXttSse5DSuc/Qe4pY8k5cTTpnY8LuSZZbEWuiX+Nlcz0xDKQlirwPaK6qgpKF/ydirX/JvWoXmT0HBr2GquRa5oL28TeJcFzLVXF69jx1lRq95aQe+p4ck8ZhySF5nibbzHNiSUYF4RbNKd1tYgnmQ4THiStU+iQyApCmebIEkyMFa4sPpBcvDu/Z/+3q8gZNJL0Ln3oePWTIb0Wm8g1zZklmBgqXFnMTa+t9pWxXPUOu957BklNJ+v4s/Ckt2qQXASYYMMh08xZgomBuwrX8NLizShQW7GLne88zv6NS0k/ZiBtL7gxZEVuTpqHovtGuBOsMVFkCcZBhSuLmfLqKvybKKI1Xra+MIXait3knXsN2QMvDNlc3oZEJpG4kmBEZCowEqgGNgK/VNXdbsTilMKVxUx5zZdctMaLJKcgySnkDfsVKflHk1rQNeQae0JkEk1S06c4YgHQx1/y4SvgdpfiiLq7CtfQ/fa3ufHVVdQpVG39mh+eu46KL94HIOu400OSi4glF5OYXOnBqOq7AYeLgcvciCPaAlfkal0tZYtfZ/fHL+PJysPTqk3I+RkpHh4YbVX9TeKKhzmYXwGvHuyHzWlfpFeWfA/Ul7GcRlXxF2QedwZtfvabkIlc2zLEtASOJZhI9kUSkTuBGuClg92nOe2LVKu+8Kp/3ED1jm9pe9FNZPUeFrIfkU3kmpbCsQTT1L5IInIlcBFwjqrGdeKIRGlpKVUbl5LWfQhZPz2N9KP7hrz9LAITTrK5FtNyuPUUaQRwK3Cmqu5zI4ZoWrhwIVdeeSWlpbtpd82zeNJbNUguNoFrWiq3niI9AWQDC0RklYj83aU4jkhlZSVTpkxh+PDh5OTk8OnHH3LlsN54/EMij4glF9OiufUUqYcb7UZTVVUVQ4cOZfXq1Vx77bU89NBDZGZmMnAgllCM8XOrB9Ns1U8XpaWlMX78eObOncsTTzxBZmamy5EZE38swRyC4uJiRowYwQcffADArbfeygUXWOEnYw7GEkyEXn/9dfr27ctHH33E1q1b3Q7HmGbBEkwTysrKmDhxImPGjKFHjx6sWrWKcePGuR2WMc2CJZgm/Otf/+LFF1/k7rvv5uOPP6Znz55uh2RMsxEPrwrEHa/Xy7p16+jXrx9XX301J510Ev3793c7LGOaHevBBFm/fj2nnHIKw4YNY9euXSQlJVlyMeYwWYLxU1WeeuopBg4cyKZNm3j66afJy8tzOyxjmjUbIuFbNHfZZZcxZ84chg8fzvPPP0/Hjh3dDsuYZs96MPgWzbVv357HHnuMefPmWXIxJkpabILZt28f119/PWvXrgVgxowZ3HDDDSQltdhfiTFR1yKHSMuXL2fChAmsX7+enj17cvzxx7sdkjEJqUX977q2tpYHHniAoUOHUl5ezqJFi/jtb3/rdljGJKwWlWCefPJJ7rjjDkaPHk1RURFnn3222yEZk9ASfoikqpSUlFBQUMCkSZPo2LEjo0aNCiljaYyJPld6MCLyJxEp8hebeldEHHlsU1payrhx4xgyZAhlZWWkpaUxevRoSy7GxIhbQ6SpqtpPVU8A5gB/iHYDixYtol+/frzxxhtMnjyZrKysaDdhjGmCKwlGVcsCDrOAqBX99nq93HzzzZx77rm0atWKxYsXc9ttt+HxeJq+2BgTVa5N8orI/SLyPTCBKPZgPB4PK1as4De/+Q0rVqxg0KBB0bq1MeYQiVM7hkSyL5L/vNuBdFW95yD3Cdx4bdB3333XZNvV1dWkpqYeVtzGmKaJyHJVHdzkeW5vSSQiXYG5qtqnqXMHDx6sy5Yti0FUxpjGRJpg3HqKFFi16WLgSzfiMMY4y611MH8WkV5AHfAdMNmlOIwxDnJrX6Sfu9GuMSa2WtSrAsaY2LIEY4xxjOtPkQ6FiOzAN2fTlHygxOFwIhUvscRLHBA/sVgcoSKNpauqFjR1UrNKMJESkWWRPEKLhXiJJV7igPiJxeIIFe1YbIhkjHGMJRhjjGMSNcFMdzuAAPESS7zEAfETi8URKqqxJOQcjDEmPiRqD8YYEwcswRhjHJOwCSZWZTkjiGOqiHzpj+VNEWntRhz+WMaIyFoRqRORmD8WFZERIrJeRDaIyG2xbj8gjmdFZLuIfO5WDP44uojIv0Vknf+/yw0uxZEuIp+JyGp/HPdF7eaqmpBfQE7A5+uBv7sUx8+AZP/nB4EHXfydHAf0At4HBse4bQ+wETgWSAVWA71d+j2cAQwEPnfrv4U/jqOAgf7P2cBXbvxOAAFa+T+nAEuAodG4d8L2YNTBspyHGMe7qlrjP1wMdHYjDn8s61R1vUvNnwhsUNVNqloN/Au4xI1AVPU/QKkbbQfFsVVVV/g/7wXWAZ1ciENVtdx/mOL/isq/l4RNMOBcWc4j8CvgHbeDcEkn4PuA4y248I8pXolIN2AAvt6DG+17RGQVsB1YoKpRiaNZJxgRWSgin4f5ugRAVe9U1S7AS8B1bsXhP+dOoMYfi2MiicUl4faKsTUSgIi0AmYBNwb1vGNGVWvVt8tHZ+BEEWmywmQkmvXGa6p6boSnvgzMBcLW/XU6DhG5ErgIOEf9A12nHMLvJNa2AF0CjjsDP7gUS9wQkRR8yeUlVX3D7XhUdbeIvA+MAI54ErxZ92AaEy9lOUVkBHArcLGq7nMjhjixFOgpIseISCowDnjL5ZhcJb4dAJ8B1qnqNBfjKKh/uikiGcC5ROnfS8Ku5BWRWfiemBwoy6mqxS7EsQFIA3b6v7VYVV0pESoio4C/AgXAbmCVqp4Xw/YvAB7D90TpWVW9P1ZtB8XxCjAMX2mCbcA9qvqMC3GcBnwIrMH39xTgDlV9O8Zx9AP+ge+/SxLwmqr+MSr3TtQEY4xxX8IOkYwx7rMEY4xxjCUYY4xjLMEYYxxjCcYY45hmvdDORJeItAUW+Q87ALXADv/xif53iGId03zgMv+7OqaZscfUJiwRuRcoV9WHg74v+P7e1IW9MHrtx6Qd4ywbIpkmiUgP//tMfwdWAF1EZHfAz8eJyAz/5/Yi8oaILPPXGBka5n5X+2vjzPfXh7nrIO0cJSJbAlaZ/tJfV2e1iDwXaXvGPTZEMpHqDfxSVSeLSGN/bx4HHlLVxf43hOcA4V6cO9H//WpgqYjMAcoD2wHwdWRARPrje+XiFFUtFZE2h9iecYElGBOpjaq6NILzzgV61ScGIE9EMlR1f9B581V1F4CIFAKnAfMaaeds4FVVLQWo//MQ2jMusARjIlUR8LmOhuUX0gM+C5FNCAdP/tUfVwSfGHDfcBOGkbZnXGBzMOaQ+Sded4lITxFJAkYF/HghcG39gYiccJDb/ExEWotIJr7Kdh830exCYFz90ChgiBRpe8YFlmDM4boV35BmEb5aL/WuBU71T8Z+AUw6yPUf4avTsxJ4RVVXNdaYqhYBDwH/8Vdem3qI7RkX2GNqE3MicjXQR1VvdDsW4yzrwRhjHGM9GGOMY6wHY4xxjCUYY4xjLMEYYxxjCcYY4xh1HHvwAAAAC0lEQVRLMMYYx/x/Rue89RItWxcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 288x216 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#还可以观察预测值与真值的散点图\n",
    "plt.figure(figsize=(4, 3))\n",
    "plt.scatter(y_train, y_train_pred_lr)\n",
    "plt.plot([-3, 3], [-3, 3], '--k')   #数据已经标准化，3倍标准差即可\n",
    "plt.axis('tight')\n",
    "plt.xlabel('True price')\n",
    "plt.ylabel('Predicted price')\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "D:\\Anaconda2\\lib\\site-packages\\sklearn\\utils\\validation.py:578: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples, ), for example using ravel().\n",
      "  y = column_or_1d(y, warn=True)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([-5.16336377e-03,  1.94769991e-04,  4.77263415e-03,  2.49912490e-03,\n",
       "       -2.29223897e-05,  4.47737676e-05,  1.30963153e-04, -1.05738010e-04,\n",
       "       -8.09577449e-05,  3.27819049e-04, -3.71207029e-05, -3.65820509e-05,\n",
       "        3.54304707e-01,  7.98616838e-01])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 线性模型，随机梯度下降优化模型参数\n",
    "# 随机梯度下降一般在大数据集上应用，其实本项目不适合用\n",
    "from sklearn.linear_model import SGDRegressor\n",
    "\n",
    "# 使用默认配置初始化线\n",
    "sgdr = SGDRegressor(max_iter=1000)\n",
    "\n",
    "# 训练：参数估计\n",
    "sgdr.fit(X_train, y_train)\n",
    "\n",
    "# 预测\n",
    "#sgdr_y_predict = sgdr.predict(X_test)\n",
    "\n",
    "sgdr.coef_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The value of default measurement of SGDRegressor on test is 0.9999998782678633\n",
      "The value of default measurement of SGDRegressor on train is 0.9999999094245067\n"
     ]
    }
   ],
   "source": [
    "# 使用SGDRegressor模型自带的评估模块(评价准则为r2_score)，并输出评估结果\n",
    "print 'The value of default measurement of SGDRegressor on test is', sgdr.score(X_test, y_test)\n",
    "print 'The value of default measurement of SGDRegressor on train is', sgdr.score(X_train, y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "作业3：岭回归（L2正则回归）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The r2 score of RidgeCV on test is 0.9999999953385482\n",
      "The r2 score of RidgeCV on train is 0.9999999968639497\n"
     ]
    }
   ],
   "source": [
    "#岭回归／L2正则\n",
    "#class sklearn.linear_model.RidgeCV(alphas=(0.1, 1.0, 10.0), fit_intercept=True, \n",
    "#                                  normalize=False, scoring=None, cv=None, gcv_mode=None, \n",
    "#                                  store_cv_values=False)\n",
    "from sklearn.linear_model import  RidgeCV\n",
    "\n",
    "#设置超参数（正则参数）范围\n",
    "alphas = [ 0.01, 0.1, 1, 10,100]\n",
    "#n_alphas = 20\n",
    "#alphas = np.logspace(-5,2,n_alphas)\n",
    "\n",
    "#生成一个RidgeCV实例\n",
    "ridge = RidgeCV(alphas=alphas, store_cv_values=True)  \n",
    "\n",
    "#模型训练\n",
    "ridge.fit(X_train, y_train)    \n",
    "\n",
    "#预测\n",
    "y_test_pred_ridge = ridge.predict(X_test)\n",
    "y_train_pred_ridge = ridge.predict(X_train)\n",
    "\n",
    "\n",
    "# 评估，使用r2_score评价模型在测试集和训练集上的性能\n",
    "print 'The r2 score of RidgeCV on test is', r2_score(y_test, y_test_pred_ridge)\n",
    "print 'The r2 score of RidgeCV on train is', r2_score(y_train, y_train_pred_ridge)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEKCAYAAAA4t9PUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xt4HfV95/H3x5Il2/IVSzbGsrEdDMRAACPbJCRpEkJi0iSGBAeDaWiXLk1bnu1eulvydEl22fbZZbttns2GXkhJS7C5xQmJkkAghFyaNJItAwabW4SxLr5g+W7Llm1J3/3jjMxBHFlHSKOjy+f1POfRnJnfzHzPsXw+mt/M+Y0iAjMzs3dqTKELMDOz4c1BYmZm/eIgMTOzfnGQmJlZvzhIzMysXxwkZmbWLw4SMzPrFweJmZn1i4PEzMz6pbjQBQyG8vLymDdvXqHLMDMbVjZu3LgnIip6azcqgmTevHnU1dUVugwzs2FFUkM+7dy1ZWZm/eIgMTOzfnGQmJlZvzhIzMysXxwkZmbWLw4SMzPrFweJmZn1i4PEzGwEer75AF97+jccbjuZ+r4cJGZmI9A//Wobf//zrYyRUt+Xg8TMbITZ13qCHz6/k88snk1ZafoDmDhIzMxGmHUbmzjR0cnqZWcPyv4cJGZmI0hnZ7C2tpEl86Zx3pmTBmWfDhIzsxHkl/V7aNh7lJsuH5yjEXCQmJmNKGtrGzijrITlF545aPt0kJiZjRC7Drbx1Eu7WVlVSWlx0aDt10FiZjZCPLi+kY7OYPXSwevWAgeJmdmI0N7RyUMbGvnguRXMnT5hUPedapBIWi7pFUn1km7PsbxU0sPJ8lpJ85L5SyU9lzw2Sbo2a51tkl5Ilvm2h2ZmwFMv7eaNQ8e5adncQd93at9UkVQE3A1cBTQDGyRVR8SLWc1uAfZHxDmSVgF3AdcDm4GqiGiXNAvYJOn7EdGerPfhiNiTVu1mZsPN2toGZk0Zx0fOnzHo+07ziGQpUB8RWyPiBPAQsKJbmxXAfcn0OuBKSYqIo1mhMQ6IFOs0MxvWtu1p5V9+s4dVS+ZSXDT4ZyzS3ONsoCnreXMyL2ebJDgOAtMBJC2TtAV4AfhCVrAE8KSkjZJuTbF+M7Nh4YH1jRSNEauWzinI/tMchCXXSGHdjyx6bBMRtcAFkt4N3Cfp8YhoA66IiB2SZgA/lvRyRPzibTvPhMytAHPnDn6foZnZYGg72cG36pr42KKZzJw8riA1pHlE0gxkx2MlsKOnNpKKgSnAvuwGEfES0ApcmDzfkfzcDTxKpgvtbSLinoioioiqioqKfr8YM7Oh6PHNO9l/9OSgjauVS5pBsgFYKGm+pBJgFVDdrU01cHMyfR3wdEREsk4xgKSzgfOAbZLKJE1K5pcBHyNzYt7MbFRaU9PI/PIy3veu6QWrIbWureSKq9uAJ4Ai4BsRsUXSnUBdRFQD9wL3S6oncySyKln9/cDtkk4CncAfRcQeSQuAR5UZX78YeCAifpTWazAzG8pe2nmIjQ37+a+//W7GjEn/viM9SXWg+oh4DHis27wvZU23AStzrHc/cH+O+VuBiwe+UjOz4WdtbQMlxWP47OLKgtbhb7abmQ1DR4638+gz2/nke2YxraykoLU4SMzMhqHvPrud1hMdgzpcfE8cJGZmw0xE5uZV7541mUvnTC10OQ4SM7Ph5pnGA7y08xA3XT6X5OKjgnKQmJkNM2trGphYWsyKS7oPFlIYDhIzs2Fkf+sJfvDCTq659CwmlqZ64W3eHCRmZsPIuo3NnGjvHBIn2bs4SMzMhonOzmBtbQNVZ0/j/DMnF7qcUxwkZmbDxL++tpdte4+y+vKhNRCtg8TMbJhYU9PAtAljufrCWYUu5S0cJGZmw8Cug238+KU3+FzVHMaNLSp0OW/hIDEzGwYe2tBIR2dww9Kh1a0FDhIzsyGvvaOTh9Y38YGF5cwrLyt0OW/jIDEzG+J+8vJudh1qG1KX/GZzkJiZDXFraho4c/I4rjx/RqFLyclBYmY2hDXsbeVffrOHVUvnUFw0ND+yh2ZVZmYGwAO1jRSNEauWDL2T7F0cJGZmQ1TbyQ4eqWvio++ewZlTxhW6nB6lGiSSlkt6RVK9pNtzLC+V9HCyvFbSvGT+UknPJY9Nkq7Nd5tmZiPFjzbvYv/Rk0P2JHuX1IJEUhFwN3A1sAi4QdKibs1uAfZHxDnAV4C7kvmbgaqIuARYDvyDpOI8t2lmNiKsqWlg3vQJXPGu8kKXclppHpEsBeojYmtEnAAeAlZ0a7MCuC+ZXgdcKUkRcTQi2pP544DowzbNzIa9l3cdoq5hPzcum8uYMYW/edXppBkks4GmrOfNybycbZLgOAhMB5C0TNIW4AXgC8nyfLZJsv6tkuok1bW0tAzAyzEzGzxraxopKR7DysvmFLqUXqUZJLkiNPJtExG1EXEBsAT4oqRxeW6TZP17IqIqIqoqKir6ULaZWWG1Hm/n0We388mLZjGtrKTQ5fQqzSBpBrKjtBLY0VMbScXAFGBfdoOIeAloBS7Mc5tmZsPad5/bzpHj7UNuuPiepBkkG4CFkuZLKgFWAdXd2lQDNyfT1wFPR0Qk6xQDSDobOA/Yluc2zcyGrYhgTU0j5585icVzpxW6nLykdsPfiGiXdBvwBFAEfCMitki6E6iLiGrgXuB+SfVkjkRWJau/H7hd0kmgE/ijiNgDkGubab0GM7PB9mzTAV7aeYi/uOZCpKF9kr1LqneOj4jHgMe6zftS1nQbsDLHevcD9+e7TTOzkWJNTQNlJUVcc2nO64iGJH+z3cxsiDhw9AQ/eH4n11w6m4mlqf6dP6AcJGZmQ8S6jc2caO9k9bKh/U327hwkZmZDQGdnsLa2kcVzp7LorMmFLqdPHCRmZkPAr7fu5fU9rUN+XK1cHCRmZkPAmpoGpk4YyycumlXoUvrMQWJmVmBvHGrjyRffYOVllYwbW1TocvrMQWJmVmAPb2iiozO4cZidZO/iIDEzK6D2jk4eXN/IBxaWM7+8rNDlvCMOEjOzAnr65d3sPNjG6mXDY1ytXBwkZmYFtLa2kZmTS/nou2cWupR3zEFiZlYgjXuP8ovftLBqyVyKi4bvx/HwrdzMbJhbu76BMRKrlg79m1edjoPEzKwAjrd38K26Zq48fwazpowvdDn94iAxMyuAH23exb7WE8Pym+zdOUjMzApgTU0DZ0+fwPvPKS90Kf3mIDEzG2Sv7DrMhm37uXHpXMaMGR43rzodB4mZ2SBbW9tASdEYVlYN75PsXVINEknLJb0iqV7S7TmWl0p6OFleK2leMv8qSRslvZD8/EjWOj9Ltvlc8piR5mswMxtIrcfb+c4z2/nERWdyRllJocsZEKndgktSEXA3cBXQDGyQVB0RL2Y1uwXYHxHnSFoF3AVcD+wBPhUROyRdSOYe7dn3nVwdEXVp1W5mlpbqTTs4crx9RJxk75LmEclSoD4itkbECeAhYEW3NiuA+5LpdcCVkhQRz0bEjmT+FmCcpNIUazUzS11EsKamgfPPnMRlZ08rdDkDJs0gmQ00ZT1v5q1HFW9pExHtwEFgerc2nwWejYjjWfP+KenWukPS8D9TZWajwnNNB9iy4xCrl81lJH10pRkkud6l6EsbSReQ6e76g6zlqyPiIuADyeN3cu5culVSnaS6lpaWPhVuZpaGtbWNTCgp4ppLu/9NPbylGSTNQPYlCZXAjp7aSCoGpgD7kueVwKPA5yPita4VImJ78vMw8ACZLrS3iYh7IqIqIqoqKioG5AWZmb1TB46e4PubdnDNpbOZNG5socsZUGkGyQZgoaT5kkqAVUB1tzbVwM3J9HXA0xERkqYCPwS+GBG/6mosqVhSeTI9FvgksDnF12BmNiDWbWzmeHvnsB4uviepBUlyzuM2MldcvQQ8EhFbJN0p6dNJs3uB6ZLqgf8IdF0ifBtwDnBHt8t8S4EnJD0PPAdsB76e1mswMxsIEcEDtY1cOncqF5w1pdDlDLjULv8FiIjHgMe6zftS1nQbsDLHen8B/EUPm71sIGs0M0vbr1/by9Y9rfz1yosLXUoq/M12M7OUraltYMr4sfz2e2YVupRUOEjMzFK0+1AbT255g5WXVTJubFGhy0mFg8TMLEUPb2iivTO4cQSeZO/iIDEzS0lHZ/Dg+kauOGc6CyomFrqc1DhIzMxS8tOXd7PjYBs3LRs542rl4iAxM0vJmtoGZkwq5aOLZha6lFQ5SMzMUtC07yg/f7WFVUvmMLZoZH/UjuxXZ2ZWIA+sb0TAqqUj9yR7FweJmdkAO97ewSMbmrjy3TM5a+r4QpeTOgeJmdkA+9HmXextPTEix9XKxUFiZjbA1tY2MveMCXxw4egYedxBYmY2gF594zDrX9/HjcvmMmbMyLl51ek4SMzMBtDamgZKisaw8rLKQpcyaBwkZmYD5OiJdr7zzHauvuhMpk8sLXQ5g8ZBYmY2QKqf28Hh4+3cdPnI/iZ7dw4SM7MBEBGsqW3g3JkTqTp7WqHLGVQOEjOzAfB880E2bz/ETZefjTQ6TrJ3cZCYmQ2ANTUNTCgp4tpLZxe6lEGXapBIWi7pFUn1km7PsbxU0sPJ8lpJ85L5V0naKOmF5OdHsta5LJlfL+mrGm3Rb2ZDzsGjJ/n+8ztYcclsJo0bW+hyBl1qQSKpCLgbuBpYBNwgaVG3ZrcA+yPiHOArwF3J/D3ApyLiIuBm4P6sdf4OuBVYmDyWp/UazMzy8e1nmmk72TlqvsneXd5BIun9kn4vma6QNL+XVZYC9RGxNSJOAA8BK7q1WQHcl0yvA66UpIh4NiJ2JPO3AOOSo5dZwOSI+HVEBPBN4Jp8X4OZ2UCLCNbWNnDJnKlcOHtKocspiLyCRNKXgT8DvpjMGgus6WW12UBT1vPmZF7ONhHRDhwEpndr81ng2Yg4nrRv7mWbXTXfKqlOUl1LS0svpZqZvTO/3rqX11paR90lv9nyPSK5Fvg00AqQHC1M6mWdXOcuoi9tJF1AprvrD/qwTZIa74mIqoioqqgYHePdmNngW1vbyJTxY/nke2YVupSCyTdITiRdSQEgqSyPdZqBOVnPK4EdPbWRVAxMAfYlzyuBR4HPR8RrWe2zxx3ItU0zs0Gx+3AbT2zexXWXVTJubFGhyymYfIPkEUn/AEyV9G+Bp4Cv97LOBmChpPmSSoBVQHW3NtVkTqYDXAc8HREhaSrwQ+CLEfGrrsYRsRM4LOny5GqtzwPfy/M1mJkNqEc2NNHeGdw4Sk+ydynOp1FE/B9JVwGHgPOAL0XEj3tZp13SbcATQBHwjYjYIulOoC4iqoF7gfsl1ZM5ElmVrH4bcA5wh6Q7knkfi4jdwB8C/wyMBx5PHmZmg6qjM3hwfRPve9d03lUxsdDlFJQyPVa9NMp0ZbVFRIek88iEyeMRcTLtAgdCVVVV1NXVFboMMxtBfvLSG9xyXx1/u3oxn7hoZJ4fkbQxIqp6a5dv19YvgFJJs8l0a/0emaMCM7NRaU1NAxWTSrlq0cxCl1Jw+QaJIuIo8Bng/0XEtWS+ZGhmNuo07TvKz15tYdWSOYwt8khTeQeJpPcCq8mcBIc8z6+YmY00D65vRMANS0f3SfYu+QbJnwC3A99JTpjPB55Orywzs6HpRHsnj9Q18ZHzZ3LW1PGFLmdIyPeo4ijQSWa8rJvIfDGw97P0ZmYjzI+27GLPkROsvtxHI13yDZK1wJ8Cm8kEipnZqLS2poE5Z4zntxZ6xIwu+QZJS0R8P9VKzMyGuN+8cZja1/fxZ8vPZ8wY38GiS75B8mVJ/wj8BDjeNTMivpNKVWZmQ9Da2kbGFomVVZW9Nx5F8g2S3wPOJzPqb1fXVgAOEjMbFY6eaOfbzzRz9YWzKJ9YWuhyhpR8g+Ti5CZTZmaj0vc37eBwW/uoHi6+J/le/luT4+6GZmajxpqaRs6dOZEl86YVupQhJ98geT/wXHL/9eeTe6Y/n2ZhZmZDxfPNB3hh+0FWLzubzMDjli3fri3fF93MRq01NQ2MH1vEtYtz3pB11Mt3GPmGtAsxMxuKDh49SfWmHVxzyWwmjxtb6HKGJI82ZmZ2Gt95tpm2k50+yX4aDhIzsx5EBGtrG7l4zlQunD2l0OUMWQ4SM7Me1GzdR/3uI6we5bfS7U2qQSJpeXKlV72k23MsL5X0cLK8VtK8ZP50ST+VdETS17qt87Nkm88ljxlpvgYzG73W1jYweVwxn3rPWYUuZUhLLUgkFQF3A1eTuQnWDTm+i3ILsD8izgG+AtyVzG8D7iAzUGQuqyPikuSxe+CrN7PRruXwcZ7YsovrLpvD+JKiQpczpKV5RLIUqI+IrRFxAngIWNGtzQrgvmR6HXClJEVEa0T8kkygmJkNukfqmjjZEdzobq1epRkks4GmrOfNybycbSKiHTgITM9j2/+UdGvdIX87yMwGWEdn8EBtI+9dMJ1zZkwsdDlDXppBkusDvvvNsPJp093qZNyvDySP38m5c+lWSXWS6lpaWnot1sysy89f3c32A8d8yW+e0gySZmBO1vNKYEdPbSQVA1OAfafbaERsT34eBh4g04WWq909EVEVEVUVFb4BjZnlb01NI+UTS7lq0cxClzIspBkkG4CFkuZLKgFWAdXd2lQDNyfT1wFPR0SPRySSiiWVJ9NjgU+SuWujmdmAaN5/lJ++sptVS+ZQUuxvSOQj37G2+iwi2iXdBjwBFAHfiIgtku4E6iKiGrgXuF9SPZkjkVVd60vaBkwGSiRdA3wMaACeSEKkCHgK+Hpar8HMRp8H1zci4AafZM9bakECEBGPAY91m/elrOk2YGUP687rYbOXDVR9ZmbZTrR38vCGJj583gxmTx1f6HKGDR+3mZklnnxxF3uOnPBJ9j5ykJiZJdbUNFA5bTwfPNcX6PSFg8TMDKjffZiarfu4Yelcisb462l94SAxMwPW1jYytkhcv2RO743tLRwkZjbqHTvRwbc3NrP8wlmUTywtdDnDjoPEzEa972/awaG2dg8X/w45SMxs1Ftb28DCGRNZNv+MQpcyLDlIzGxUe6H5IJuaD7J62Vw8Buw74yAxs1FtTU0D48cWce3iykKXMmw5SMxs1Dp47CTVm3bw6YvPYsr4sYUuZ9hykJjZqPXoM80cO9nhb7L3k4PEzEaliGBNbSPvqZzCRZVTCl3OsOYgMbNRaf3r+6jffYSblvlopL8cJGY2Kq2pbWTSuGI+dfFZhS5l2HOQmNmo03L4OD/avJPPLq5kfElRocsZ9hwkZjbqfGtjEyc7gpsu9zfZB4KDxMxGlY7O4IHaRi5fcAbnzJhU6HJGhFSDRNJySa9Iqpd0e47lpZIeTpbXSpqXzJ8u6aeSjkj6Wrd1LpP0QrLOV+WvoppZH/zi1Raa9x9jtU+yD5jUgkRSEXA3cDWwCLhB0qJuzW4B9kfEOcBXgLuS+W3AHcCf5tj03wG3AguTx/KBr97MRqq1tQ2UTyzl4xecWehSRow0j0iWAvURsTUiTgAPASu6tVkB3JdMrwOulKSIaI2IX5IJlFMkzQImR8SvIyKAbwLXpPgazGwE2X7gGE+/vJvrl1RSUuye/YGS5js5G2jKet6czMvZJiLagYPA9F622dzLNs3McnqwtpEAVi3xSfaBlGaQ5Dp3Ee+gzTtqL+lWSXWS6lpaWk6zSTMbDU52dPLQhiY+fN4M5pwxodDljChpBkkzkH3PykpgR09tJBUDU4B9vWwze4jOXNsEICLuiYiqiKiqqKjoY+lmNtI8ueUN9hw57kt+U5BmkGwAFkqaL6kEWAVUd2tTDdycTF8HPJ2c+8gpInYChyVdnlyt9XngewNfupmNNGtqGpg9dTy/de6MQpcy4hSnteGIaJd0G/AEUAR8IyK2SLoTqIuIauBe4H5J9WSORFZ1rS9pGzAZKJF0DfCxiHgR+EPgn4HxwOPJw8ysR/W7j/DrrXv5zx8/j6Ix/sbAQEstSAAi4jHgsW7zvpQ13Qas7GHdeT3MrwMuHLgqzWyke6C2keIx4nNVc3pvbH3m69/MbEQ7dqKDdRub+PiFZ1IxqbTQ5YxIDhIzG9F+8PwODrW1e7j4FDlIzGxEW1PbyLsqyrh8wRmFLmXEcpCY2Yi1eftBNjUdYPWys/GwfOlxkJjZiLW2toFxY8fw2csqe29s75iDxMxGpENtJ/nuszv49MVnMWX82EKXM6I5SMxsRHr0me0cO9nh4eIHgYPEzEaciGBtbQMXzZ7CxXOmFrqcEc9BYmYjzoZt+3n1jSMeV2uQOEjMbMRZU9PApHHFfOriswpdyqjgIDGzEWXPkeM8vnknn11cyYSSVEeBsoSDxMxGlG/VNXOyI1i9zN1ag8VBYmYjRmdn8MD6BpbOP4OFMycVupxRw0FiZiPGL37TQtO+Y9x0uS/5HUwOEjMbMdbUNDK9rITlF5xZ6FJGFQeJmY0I2w8c4+mX3+BzS+ZQUuyPtsHkd9vMRoSH1zcSwI1LfZJ9sDlIzGzYO9nRyUMbmvjQuRXMOWNCocsZdVINEknLJb0iqV7S7TmWl0p6OFleK2le1rIvJvNfkfTxrPnbJL0g6TlJdWnWb2bDw49ffIPdh497XK0CSe3bOpKKgLuBq4BmYIOk6oh4MavZLcD+iDhH0irgLuB6SYuAVcAFwFnAU5LOjYiOZL0PR8SetGo3s+FlTU0Ds6eO58Pnzyh0KaNSmkckS4H6iNgaESeAh4AV3dqsAO5LptcBVypz95kVwEMRcTwiXgfqk+2Zmb3Fay1H+NfX9nLD0jkUjfHNqwohzSCZDTRlPW9O5uVsExHtwEFgei/rBvCkpI2Sbu1p55JulVQnqa6lpaVfL8TMhq4HahspHiM+t2ROoUsZtdIMklx/GkSebU637hURsRi4GvhjSR/MtfOIuCciqiKiqqKiIt+azWwYaTvZwbqNzXz8gjOZMWlcocsZtdIMkmYg+0+ESmBHT20kFQNTgH2nWzciun7uBh7FXV5mo9YPnt/JwWMnWe3h4gsqzSDZACyUNF9SCZmT59Xd2lQDNyfT1wFPR0Qk81clV3XNBxYC6yWVSZoEIKkM+BiwOcXXYGZD2JqaBhZUlPHeBdMLXcqoltpVWxHRLuk24AmgCPhGRGyRdCdQFxHVwL3A/ZLqyRyJrErW3SLpEeBFoB3444jokDQTeDRzPp5i4IGI+FFar8HMhq7N2w/yXNMB7vjkIpLPBCuQVAfrj4jHgMe6zftS1nQbsLKHdf8S+Mtu87YCFw98pWY23KytbaS0eAzXLa4sdCmjnr/ZbmbDzuG2k3zvue186uKzmDJhbKHLGfUcJGY27Dz67HaOnujwcPFDhIPEzIaViGBtTSMXzp7MxZVTCl2O4SAxs2GmrmE/r7xxmNXLzvZJ9iHCQWJmw8qamgYmlRaz4pKzCl2KJRwkZjZs7D1ynMdf2MVnFs9mQkmqF51aHzhIzGzY+NbGZk50dLLaJ9mHFAeJmQ0LnZ3BA7WNLJ13BufOnFTociyLjw3NbEg73t5B076j/OyVFhr3HeU/fezcQpdk3ThIzKzgOjuDXYfaeH1PK1tbjrB1T2sy3Urz/qN0JmN/nz19AssvPLOwxdrbOEjMbNAcPHqSrXuOnAqJ1/e0snVPK9v2tHLsZMepduPHFjG/vIz3VE7hmkvOYn5FGQvKJ3LuzEmUFhcV8BVYLg4SMxtQx9s7aNh79M2gaMkEx+t7WtnbeuJUu6IxYs608cwvL+N975rO/PIyFpSXsaBiIjMnl/o7IsOIg8TM+qyzM9h5qO1USGxtaU26o46wff+xU11RABWTSplfXsZVi2ayoKKM+eUTmV9extwzJlBS7Ot9RgIHiZn16MDRE2w91Q311i6p4+2dp9qVlRQxv6KMS+ZM49pLK3lXRRnzy8uYV17G5HEeVHGkc5CYjXJtJzNdUa/vOcJrSUh0dUntP3ryVLuiMWLuGRNYUF7G+88pP3XeYkFFGTMmuStqNHOQmI0CnZ3BjoPHTh1NvL6nldeSbqntB44RWV1RM5KuqOUXzmJBeebIYkFFGXPOmMDYIndF2ds5SMxGkP2tXV1RR7KOLFrZtvftXVELKiayeO40rrusMjnRPZH5FWVMLPXHgvVNqr8xkpYD/5fMrXb/MSL+V7flpcA3gcuAvcD1EbEtWfZF4BagA/h3EfFEPts0G+naTnawbW8rrycnuLvOX2zd08qBrK6o4jFi7vRMV9QHzy1nQcXEU1dGVbgrygZQakEiqQi4G7gKaAY2SKqOiBezmt0C7I+IcyStAu4Crpe0iMz92y8AzgKektT1ddbetmk27HV0BjsOHMtcCdXtC3o7Dr61K2rm5ExX1CcumpVcPpu5Mqpy2nh3RdmgSPOIZClQn9xnHUkPASuA7A/9FcB/S6bXAV9T5s+kFcBDEXEceF1SfbI98tim2ZAVSQJ0BcGBYye7fZM70yW1be9RTmR1RU0sLWZBRRlL5k1jfvmc5ER35qood0VZoaX5GzgbaMp63gws66lNRLRLOghMT+bXdFt3djLd2zYHzO/ft4Fte4+mtfkRKSI49cdyvPnj1AfoqXZdy+LN6ay/srtvK2d7uq+Xa9mb+31zP1m15NOet374v3W9N5fl2kau19WTsUWZq6Lml0/kQ+fNyDrRPZHyiSXuirIhK80gyfVb3/2/VE9tepqf6zg9539TSbcCtwLMnTu35ypP4+zpZR6O4Z3Qm/+AXR9+Aro+B3Mt49SyZF5WW+VY9uZnan7t31z21l8tKXe7N/ejt8zLrjNn+x62hfSWdl3bKCst4l3JuYvKaeMpdleUDUNpBkkzMCfreSWwo4c2zZKKgSnAvl7W7W2bAETEPcA9AFVVVXn+TfhWd3xy0TtZzcxsVEnzz58NwEJJ8yWVkDl5Xt2tTTVwczJ9HfB0ZPoFqoFVkkolzQcWAuvz3KaZmQ2i1I5IknMetwFPkLlU9xsRsUXSnUBdRFQD9wL3JyfT95EJBpJ2j5A5id4O/HFEdADk2mZar8HMzHqnyPdM4DBWVVUVdXV1hS7DzGxYkbQxIqp6a+cze2Zm1i8OEjMz6xcHiZmZ9YuDxMzM+sVBYmZm/TIqrtqS1AJu408AAAAHc0lEQVQ0vMPVy4E9A1jOQHFdfeO6+sZ19c1IrevsiKjordGoCJL+kFSXz+Vvg8119Y3r6hvX1TejvS53bZmZWb84SMzMrF8cJL27p9AF9MB19Y3r6hvX1Tejui6fIzEzs37xEYmZmfWLg6QbSX8l6WVJz0t6VNLUHtotl/SKpHpJtw9CXSslbZHUKanHqzAkbZP0gqTnJKU+UmUf6hrs9+sMST+W9Jvk57Qe2nUk79VzklK7JUFvrz+5ZcLDyfJaSfPSqqWPdf2upJas9+j3B6Gmb0jaLWlzD8sl6atJzc9LWpx2TXnW9SFJB7Peqy8NUl1zJP1U0kvJ/8U/ydEm3fcsIvzIegAfA4qT6buAu3K0KQJeAxYAJcAmYFHKdb0bOA/4GVB1mnbbgPJBfL96ratA79f/Bm5Ppm/P9e+YLDsyCO9Rr68f+CPg75PpVcDDQ6Su3wW+Nli/T8k+PwgsBjb3sPwTwONkbkB5OVA7ROr6EPCDwXyvkv3OAhYn05OAV3P8O6b6nvmIpJuIeDIi2pOnNWTuwtjdUqA+IrZGxAngIWBFynW9FBGvpLmPdyLPugb9/Uq2f18yfR9wTcr7O518Xn92veuAK5X+TdoL8e/Sq4j4BZn7E/VkBfDNyKgBpkqaNQTqKoiI2BkRzyTTh4GXgNndmqX6njlITu/fkEnx7mYDTVnPm3n7P1yhBPCkpI3JfeuHgkK8XzMjYidk/qMBM3poN05SnaQaSWmFTT6v/1Sb5A+Zg8D0lOrpS10An026Q9ZJmpNj+WAbyv//3itpk6THJV0w2DtPukQvBWq7LUr1PUvznu1DlqSngDNzLPrziPhe0ubPydydcW2uTeSY1+/L3/KpKw9XRMQOSTOAH0t6OflLqpB1Dfr71YfNzE3erwXA05JeiIjX+ltbN/m8/lTeo17ks8/vAw9GxHFJXyBz1PSRlOvqTSHeq3w8Q2ZIkSOSPgF8l8xtwgeFpInAt4F/HxGHui/OscqAvWejMkgi4qOnWy7pZuCTwJWRdDB20wxk/2VWCexIu648t7Ej+blb0qNkui/6FSQDUNegv1+S3pA0KyJ2Jofwu3vYRtf7tVXSz8j8NTfQQZLP6+9q0yypGJhC+t0ovdYVEXuznn6dzHnDQkvl96m/sj+8I+IxSX8rqTwiUh+DS9JYMiGyNiK+k6NJqu+Zu7a6kbQc+DPg0xFxtIdmG4CFkuZLKiFzcjS1K37yJalM0qSuaTIXDuS8wmSQFeL9qgZuTqZvBt525CRpmqTSZLocuAJ4MYVa8nn92fVeBzzdwx8xg1pXt370T5Ppfy+0auDzyZVIlwMHu7oxC0nSmV3ntSQtJfP5uvf0aw3IfgXcC7wUEX/TQ7N037PBvsJgqD+AejJ9ic8lj64rac4CHstq9wkyV0e8RqaLJ+26riXzV8Vx4A3gie51kbn6ZlPy2DJU6irQ+zUd+Anwm+TnGcn8KuAfk+n3AS8k79cLwC0p1vO21w/cSeYPFoBxwLeS37/1wIK036M86/qfye/SJuCnwPmDUNODwE7gZPK7dQvwBeALyXIBdyc1v8BprmIc5Lpuy3qvaoD3DVJd7yfTTfV81ufWJwbzPfM3283MrF/ctWVmZv3iIDEzs35xkJiZWb84SMzMrF8cJGZm1i8OErPTkHSkn+uvS741f7o2P9NpRk7Ot0239hWSfpRve7P+cJCYpSQZa6koIrYO9r4jogXYKemKwd63jT4OErM8JN8I/itJm5W538v1yfwxyVAYWyT9QNJjkq5LVltN1jfqJf1dMkDkFkn/vYf9HJH015KekfQTSRVZi1dKWi/pVUkfSNrPk/QvSftnJL0vq/13kxrMUuUgMcvPZ4BLgIuBjwJ/lQwf8hlgHnAR8PvAe7PWuQLYmPX8zyOiCngP8FuS3pNjP2XAMxGxGPg58OWsZcURsRT491nzdwNXJe2vB76a1b4O+EDfX6pZ34zKQRvN3oH3kxkFtwN4Q9LPgSXJ/G9FRCewS9JPs9aZBbRkPf9cMrR/cbJsEZlhLbJ1Ag8n02uA7AH4uqY3kgkvgLHA1yRdAnQA52a1301mqBqzVDlIzPLT002mTnfzqWNkxtBC0nzgT4ElEbFf0j93LetF9hhGx5OfHbz5f/c/kBnj7GIyPQxtWe3HJTWYpcpdW2b5+QVwvaSi5LzFB8kMrvhLMjd+GiNpJpnbrXZ5CTgnmZ4MtAIHk3ZX97CfMWRG/wW4Mdn+6UwBdiZHRL9D5va5Xc5laIz+bCOcj0jM8vMomfMfm8gcJfyXiNgl6dvAlWQ+sF8lc2e6g8k6PyQTLE9FxCZJz5IZHXYr8Kse9tMKXCBpY7Kd63up62+Bb0taSWZ03tasZR9OajBLlUf/NesnSRMjc1e86WSOUq5IQmY8mQ/3K5JzK/ls60hETBygun4BrIiI/QOxPbOe+IjErP9+IGkqUAL8j4jYBRARxyR9mcy9sRsHs6Ck++1vHCI2GHxEYmZm/eKT7WZm1i8OEjMz6xcHiZmZ9YuDxMzM+sVBYmZm/eIgMTOzfvn/7R7D6SwexOkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('alpha is:', 0.01)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>coef_lr</th>\n",
       "      <th>coef_ridge</th>\n",
       "      <th>columns</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>[0.7992015607529477]</td>\n",
       "      <td>[0.7990593077353303]</td>\n",
       "      <td>registered</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>[0.35431900383432696]</td>\n",
       "      <td>[0.35431653529558116]</td>\n",
       "      <td>casual</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>[2.664146216407351e-15]</td>\n",
       "      <td>[-0.00015264733927639051]</td>\n",
       "      <td>instant</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>[4.306471171686414e-16]</td>\n",
       "      <td>[2.928478783534305e-05]</td>\n",
       "      <td>temp</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>[9.683329465026209e-17]</td>\n",
       "      <td>[-1.1450215976205325e-05]</td>\n",
       "      <td>hum</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>[6.768049791686516e-18]</td>\n",
       "      <td>[9.790727595812565e-06]</td>\n",
       "      <td>weekday</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>[-1.3975909410059654e-17]</td>\n",
       "      <td>[-2.4417945238636163e-05]</td>\n",
       "      <td>weathersit</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>[-8.086278779127036e-17]</td>\n",
       "      <td>[-1.1967389549460278e-05]</td>\n",
       "      <td>windspeed</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>[-1.9648465161928974e-16]</td>\n",
       "      <td>[-3.5643288828907105e-06]</td>\n",
       "      <td>holiday</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>[-2.6278332966459976e-16]</td>\n",
       "      <td>[5.1746401972058154e-05]</td>\n",
       "      <td>season</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>[-3.1207395989199734e-16]</td>\n",
       "      <td>[3.367550173681566e-05]</td>\n",
       "      <td>workingday</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>[-3.7441501753193594e-16]</td>\n",
       "      <td>[2.770302526397539e-05]</td>\n",
       "      <td>atemp</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>[-7.378407600702821e-16]</td>\n",
       "      <td>[6.409393514759942e-05]</td>\n",
       "      <td>mnth</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>[-2.2732219189246263e-15]</td>\n",
       "      <td>[0.0002087671911854272]</td>\n",
       "      <td>yr</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                      coef_lr                 coef_ridge     columns\n",
       "13       [0.7992015607529477]       [0.7990593077353303]  registered\n",
       "12      [0.35431900383432696]      [0.35431653529558116]      casual\n",
       "0     [2.664146216407351e-15]  [-0.00015264733927639051]     instant\n",
       "8     [4.306471171686414e-16]    [2.928478783534305e-05]        temp\n",
       "10    [9.683329465026209e-17]  [-1.1450215976205325e-05]         hum\n",
       "5     [6.768049791686516e-18]    [9.790727595812565e-06]     weekday\n",
       "7   [-1.3975909410059654e-17]  [-2.4417945238636163e-05]  weathersit\n",
       "11   [-8.086278779127036e-17]  [-1.1967389549460278e-05]   windspeed\n",
       "4   [-1.9648465161928974e-16]  [-3.5643288828907105e-06]     holiday\n",
       "1   [-2.6278332966459976e-16]   [5.1746401972058154e-05]      season\n",
       "6   [-3.1207395989199734e-16]    [3.367550173681566e-05]  workingday\n",
       "9   [-3.7441501753193594e-16]    [2.770302526397539e-05]       atemp\n",
       "3    [-7.378407600702821e-16]    [6.409393514759942e-05]        mnth\n",
       "2   [-2.2732219189246263e-15]    [0.0002087671911854272]          yr"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mse_mean = np.mean(ridge.cv_values_, axis = 0)\n",
    "plt.plot(np.log10(alphas), mse_mean.reshape(len(alphas),1)) \n",
    "\n",
    "#这是为了标出最佳参数的位置，不是必须\n",
    "#plt.plot(np.log10(ridge.alpha_)*np.ones(3), [0.28, 0.29, 0.30])\n",
    "\n",
    "plt.xlabel('log(alpha)')\n",
    "plt.ylabel('mse')\n",
    "plt.show()\n",
    "\n",
    "print ('alpha is:', ridge.alpha_)\n",
    "\n",
    "# 看看各特征的权重系数，系数的绝对值大小可视为该特征的重要性\n",
    "fs = pd.DataFrame({\"columns\":list(columns), \"coef_lr\":list((lr.coef_.T)), \"coef_ridge\":list((ridge.coef_.T))})\n",
    "fs.sort_values(by=['coef_lr'],ascending=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "L1正则：lasso回归"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#### Lasso／L1正则\n",
    "# class sklearn.linear_model.LassoCV(eps=0.001, n_alphas=100, alphas=None, fit_intercept=True, \n",
    "#                                    normalize=False, precompute=’auto’, max_iter=1000, \n",
    "#                                    tol=0.0001, copy_X=True, cv=None, verbose=False, n_jobs=1,\n",
    "#                                    positive=False, random_state=None, selection=’cyclic’)\n",
    "from sklearn.linear_model import LassoCV\n",
    "\n",
    "#设置超参数搜索范围\n",
    "#alphas = [ 0.01, 0.1, 1, 10,100]\n",
    "\n",
    "#生成一个LassoCV实例\n",
    "#lasso = LassoCV(alphas=alphas)  \n",
    "lasso = LassoCV()  \n",
    "\n",
    "#训练（内含CV）\n",
    "lasso.fit(X_train, y_train)  \n",
    "\n",
    "#测试\n",
    "y_test_pred_lasso = lasso.predict(X_test)\n",
    "y_train_pred_lasso = lasso.predict(X_train)\n",
    "\n",
    "\n",
    "# 评估，使用r2_score评价模型在测试集和训练集上的性能\n",
    "print 'The r2 score of LassoCV on test is', r2_score(y_test, y_test_pred_lasso)\n",
    "print 'The r2 score of LassoCV on train is', r2_score(y_train, y_train_pred_lasso)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
